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ABSTRACT 

Rotating neutron stars, or pulsars, are plausibly the source of power behind many astrophysi¬ 
cal systems, such as gamma-ray bursts, supernovae, pulsar wind nebulae and supernova rem¬ 
nants. In the past several years, 3D numerical simulations made it possible to compute pulsar 
spindown luminosity from first principles and revealed that oblique pulsar winds are more 
powerful than aligned ones. However, what causes this enhanced power output of oblique pul¬ 
sars is not understood. In this work, using time-dependent 3D magnetohydrodynamic (MHD) 
and force-free simulations, we show that, contrary to the standard paradigm, the open mag¬ 
netic flux, which carries the energy away from the pulsar, is laterally non-uniform. We argue 
that this non-uniformity is the primary reason for the increased luminosity of oblique pul¬ 
sars. To demonstrate this, we construct simple analytic descriptions of aligned and orthogonal 
pulsar winds and combine them to obtain an accurate 3D description of the pulsar wind for 
any obliquity. Our approach describes both the warped magnetospheric current sheet and the 
smooth variation of pulsar wind properties outside of it. We find that generically the magne¬ 
tospheric current sheet separates plasmas that move at mildly relativistic velocities relative 
to each other. This suggests that the magnetospheric reconnection is a type of driven, rather 
than free, reconnection. The jump in magnetic field components across the current sheet de¬ 
creases with increasing obliquity, which could be a mechanism that reduces dissipation in 
near-orthogonal pulsars. Our analytical description of the pulsar wind can be used for con¬ 
structing models of pulsar gamma-ray emission, pulsar wind nebulae, and magnetar-powered 
core-collapse gamma-ray bursts and supernovae. 
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1 INTRODUCTION 

A rotating magnet in vacuum is the simplest description of a ro¬ 
tating neutron star, or pulsar. The magnet loses its energy to mag¬ 
netodipole radiation (Deutsch 1955). However, real pulsar magne¬ 
tospheres are likely filled with plasma (Goldreich & Julian 1969), 
which supports charges and currents that modify the pulsar spin- 
down. For a long time it has been realized that plasma-filled mag¬ 
netospheres open up to infinity due to relativistic rotation (Ingra¬ 
ham 1973; Michel 1974). The open magnetic flux carries away the 
spindown luminosity of the star in the form of a pulsar wind. 

Whereas the origin of spindown luminosity of vacuum rota¬ 
tors has been well-understood due to the availability of the analytic 
model, the spindown of plasma-filled magnetospheres can only be 
studied numerically, and progress became possible only recently. 
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The primary difficulty in computing pulsar magnetosphere struc¬ 
ture is in connecting the conditions on the star to the properties of 
the emerging pulsar wind. Of particular importance is the determi¬ 
nation of the twist of magnetic field lines at the stellar surface, or 
the magnitude of the electric current flowing through the magneto¬ 
sphere. Early pulsar magnetosphere solutions prescribed the mag¬ 
nitude of the twist, and. as a result, the solutions extended smoothly 
only to the light cylinder radius (Beskin, Gurevich & Istomin 1993; 
Mestel & Shibata 1994; Beskin 1997), the distance from the rota¬ 
tional axis at which a field line rigidly rotating at the frequency of 
the central object would move at the speed of light. 


where £2 = In/P is pulsar angular frequency and P is the period. 

In real magnetospheres, the magnitude of the magnetospheric 
current is set self-consistently by the large-scale properties of the 
magnetosphere. The rotation of the star twists the magnetic field 
lines, which sends out a torsional Alfven wave that carries the twist 
outward. On the way out, the Alfven wave undergoes partial re- 
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flections that reach back to the star and modify the surface value 
of the magnetic field twist, which determines the strength of the 
magnetospheric current. It is accounting for this communication 
by Alfven waves and determining the correct value of the mag¬ 
netospheric current that was the challenge. Contopoulos, Kazanas 
& Fendt (1999) found an iterative way to determine the strength 
of the magnetospheric current self-consistently. They obtained the 
first solution for an aligned force-free pulsar magnetosphere that 
extended out to infinity. Their results were subsequently verified by 
other groups within force-free and magnetohydrodynamic (MHD) 
approximations (e.g., Gmzinov 2005; Timokhin 2006; McKin¬ 
ney 2006; Komissarov 2006; Parfrey, Beloborodov & Hui 2012; 
Ruiz, Paschalidis & Shapiro 2014) as well as using a particle-in- 
cell (PIC) approach (Philippov & Spitkovsky 2014; Chen & Be¬ 
loborodov 2014; Cerutti, Philippov, Parfrey & Spitkovsky 2014; 
Belyaev 2014). 

Spitkovsky (2006) carried out the first 3D, oblique pul¬ 
sar magnetosphere simulations. Using the force-free approxima¬ 
tion, he found that pulsar spindown luminosity increases with in¬ 
creasing obliquity angle, a , the angle between the rotational and 
magnetic axes. More recently, these results were confirmed us¬ 
ing time-dependent 3D force-free (Kalapotharakos & Contopoulos 
2009; Petri 2012; Kalapotharakos, Contopoulos & Kazanas 2012), 
MHD (Tchekhovskoy, Spitkovsky & Li 2013), and PIC (Philippov, 
Spitkovsky & Cerutti 2014) studies. 

Despite these advances, what sets the spindown luminosity of 
plasma-filled pulsars remains elusive. Understanding this and con¬ 
structing a detailed, 3D model of a pulsar wind are the goals of this 
work. We start with a discussion of pulsar spindown in the split- 
monopole approximation and confront it with the numerical simu¬ 
lations in §2, discuss the effects of pulsar wind non-uniformity on 
pulsar spindown rate and develop a toy model of a non-uniform pul¬ 
sar wind in §3. We construct an analytic description of aligned and 
orthogonal plasma-filled pulsar winds in §4. We then try and get an 
inspiration on how to construct solutions for intermediate inclina¬ 
tion angles in §5 and construct our analytical description of the pul¬ 
sar wind for an arbitrary inclination angle in §6. We discuss our re¬ 
sults and their observational implications in §7. We will use spheri¬ 
cal polar coordinates, (r, 6 , </>), cylindrical coordinates, (R , z, <p), and 
Cartesian coordinates, (x,y, z). 


2 DEPENDENCE OF SPIN-DOWN POWER ON 

INCLINATION AND THE STRUCTURE OF PULSAR 
WIND 

Our goal is to understand the physics of spindown of plasma-filled 
pulsars. The solid blue line in Fig. 1 shows the dependence of 
pulsar spindown on inclination (see Spitkovsky 2006; Petri 2012; 
Tchekhovskoy et al. 2013; Philippov et al. 2014). Clearly, oblique 
pulsars have higher spindown rates, but the reason for this increase 
is unclear. In the widely used and highly successful standard model 
of the asymptotic pulsar wind, the wind structure is radial and is 
that of a split-monopole outflow (Bogovalov 1999). The wind in¬ 
deed appears to be radial far from the pulsar, as is seen in Fig. 2. 
However, this brings up a problem: the spindown of split-monopole 
solution is independent of inclination! In it, \r 2 B r \ = constant and 
is time-independent. In fact, as we will derive below (see eq. 14), 
the spindown luminosity is set by the rotational frequency of the 
star, fl* = In/P, and the amount of open magnetic flux in the split- 



Figure 1. Blue solid line shows the dependence of spindown luminosity 
on the inclination angle, a, in a series of MHD simulations (Tchekhovskoy 
et al. 2013). Red dotted line shows the corresponding spindown luminosity 
for a split-monopole wind (same open magnetic flux as in the simulation). It 
accounts for about 40% of the increase but falls short of explaining the full 
magnitude of the increase in luminosity with increasing a. Green dashed 
line shows the spindown for a non-uniform distribution of open magnetic 
flux (same non-uniformity as in the simulation). It matches the simulated lu¬ 
minosity well. The non-uniformity of the open magnetic flux plays a crucial 
role in setting the spindown luminosity of oblique pulsars. 



X [Rlc] 


Figure 2. Magnetic field structure of pulsar wind in the fl-Q. plane of an 
oblique, a = 60°, MHD simulation. Solid black lines show magnetic field 
lines and color shows rB y , into-the-plane magnetic field component (which 
we rescaled by the r prefactor to remove the trivial radial trend). That the 
magnetic field lines are predominantly straight suggests that the magneto¬ 
spheric outflow is nearly radial. However, as indicated by non-trivial varia¬ 
tions of rBy, whose absolute magnitude would be a constant along 9 = con¬ 
stant in a split-monopole solution, this does not mean that the solution is 
well-described by a split-monopole (see also Fig. 4). 


monopole solution, 

®open = d/2) fj\B,\d0J, (2) 

where the integral is over both hemispheres and the factor of (1/2) 
converts it to a single hemisphere. For a split-monopole magnetic 
field, eq. (2) reduces to ® 0 p en = 27ir 2 |B,.|, and the spindown lumi- 
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Figure 3. The square of the open magnetic flux <D 0 pen measured according 
to eq. (2) at a sphere of radius r = 2R^c, plotted versus the inclination 
angle, a, in a series of MHD simulations. The open magnetic flux increases 
with increasing inclination angle. This increase accounts for « 40% of the 
enhancement in spindown luminosity of oblique rotators. The remaining 
60% come from the redistribution of the open magnetic flux toward the 
equatorial plane (see §3). 


nosity is 

£ 2 2 ® 2 

L = ° pen , (uniform \rB r \) (3) 

6 7T~C 

where c is the speed of light. 

A natural possibility is that the amount of open magnetic flux 
depends on the inclination angle, ® open = ® 0 pen(aT an d this causes 
the dependence of the spindown luminosity on a. Fig. 3 shows that 
while this is indeed the case to a degree, the dependence fl> 0 pen(a) is 
too weak to fully account for the enhancement in spindown rate for 
oblique pulsars. In fact, it accounts for about 40% of the luminos¬ 
ity enhancement, as shown with red dotted lines in Figs. 1 and 3. 
What is the cause of the rest of the enhancement, i.e., the difference 
between the dotted red and solid blue lines? 


3 NON-UNIFORMITY OF ASYMPTOTIC PULSAR WIND 
AND SPINDOWN RATE 

As the colour variation in Fig. 4 suggests, the product |r 2 B,.| is not 
constant in the simulated pulsar wind. However, the split-monopole 
model of the pulsar wind and the associated spindown luminosity 
given by equation (3) assume that |r 2 B,-| is a constant everywhere, 
i.e., apart from B, changing sign at current sheet crossings, B, is 
laterally-uniform (i.e., it has no dependence on 9 and i p) and follows 
a simple power-law dependence in radius. According to Fig. 4, this 
approximation is at best very rough: B r shows variations away from 
the expected scaling, 1 /r 2 , which indicates the presence of a sub¬ 
stantial, order unity non-uniformity in the pulsar wind. Can these 
non-uniformities lead to deviations from the spindown rate given 
by equation (3)? In section §3.1 we will demonstrate that this is 
indeed the case. 


3.1 Velocity and freeze-in condition in ideal MHD 

To describe the asymptotic structure of pulsar wind, we will make 
use of the fact that the wind structure is accurately described by 
a 2D distribution of just a single variable, B r , in 8-ip plane. This 



Figure 4. The same magnetospheric structure as Fig. 2 but with colour 
showing scaled radial magnetic field, r 1 B, . The r 2 -prefactor eliminates 
the obvious radial trend: for a perfectly monopolar field, we would have 
r 2 B r = const. The substantial deviations from constancy of r 2 B, are a 
strong indication that the radial magnetic field in the pulsar magnetosphere 
is substantially non-uniform and deviates from the simple monopolar de¬ 
pendence, B r oc 1/r 2 . 


can be seen as follows. In ideal MHD/force-free a pulsar magne¬ 
tosphere rotates at the stellar angular frequency. Thus, the velocity 
can be decomposed into motion along magnetic field lines with ve¬ 
locity t)|| and together with the magnetic field lines corotating with 
the star at v v = Qr sin 9: 

„ 2 

IT = U|i -57 + Dr sin (4) 

N 

Due to the ideal MHD condition, the electric field is 2 = —(i f/c) x 
B = -{Q.rs'm9/c)B r <? g . 

Asymptotically far from the pulsar, rsind » R LC , the mag¬ 
netic field is toroidally-dominated, B v » B ,, B g and moves outward 
at a drift velocity, 

_ c\2x2\ _ c E e B „ 
dr B 2 ~ B 2 ' 

Since at such large distances the wind moves relativistically fast, 
"Udr ~ c, we have: 


B^ x -Eg = -(fir sin 8/c)B r . ( 6 ) 

Using eq. ( 6 ), for a given distribution of B r , we can reconstruct the 
rest of the field components. 


3.2 Poynting flux 


Using eq. ( 6 ), we can write the Poynting flux, which is carried by 
the outflowing toroidal magnetic field, as 

S r = ^-(2 x B) ~ ^—B 2 r 2 sin 2 9. (7) 

4n 4/re 

Thus, the spindown luminosity is given by 


■ JJ " S r du) = 2n (S r ) v r 2 sin 8d8 » — J' 


(B 2 ) v r 4 sin 3 9d9, 
( 8 ) 
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where is the averaging over <p, and we used the expression for 
the ip- average Poynting flux, 

I dL\ Or , , , 

<Sr> Had = 4 Vc {B ^ r * m ~ e ' (9) 

that trivially follows from the definition of S r , eq. (7). Both L and 
its angular distribution < dL/da ;) v are fully determined by the sur¬ 
face distribution of a single quantity, B r (8,<p), or, more precisely, 
by the </)-average of its square, (B 2 ) v . 

If the magnetic field is split-monopolar, i.e., |B,.| = BoUo/f ) 2 
and Bo = © op en/2 nr^, eq. (8) reduces to the standard split- 
monopolar expression, eq. (3). Note, however, that the real pulsar 
wind is non-uniform laterally, and this approximation does not cap¬ 
ture this nonuniformity. However, even despite the non-uniformity, 
the pulsar wind remains very nearly radial. This means that the ra¬ 
dial scaling of B, is a trivial power-law, oc r~ 2 . 


3.3 Dipolar approximation for radial magnetic field 

To make progress, we will approximate the angular distribution of 
B, in the open magnetic flux region by that of a dipole ft inclined 
at an angle a in the x—z, or tp = 0, plane. As we will see below, 
this is an excellent approximation for near-orthogonal rotators and 
a reasonable approximation for other obliquities. We obtain: 

B, = B 0 (r/r 0 y 2 cos 8 m , (10) 

where 9 m is magnetic colatitude, the angle away from the magnetic 
axis, 

6 m = arccosfsin tr sin 0 cos ^ + cos # cos or), ( 11 ) 

and 

Bq = ^open/^. (12) 

The angular distribution of B, due to eq. (10) is clearly nonuniform: 
B r is strongest at the magnetic poles (cos0 m = ±1) and weakest 
along the magnetic equator (cos 8 m = 0 ). 

To find out how this nonuniformity affects the spindown, we 
first compute the angular distribution of </>-averaged B 2 : 

{B 2 ) v = Bg(r/r 0 y 4 ( 0.5 sin 2 a sin 2 8 + cos 2 0 cos 2 a). (13) 

Now, by plugging eq. (13) into eq. ( 8 ) and integrating, we obtain : 1 

Q 2 (I> 2 

AtipoieCff) = + sin 2 a). (14) 

This is very close to the cr-dependence that is observed in numerical 
simulations and shown in Fig. 1 with the solid blue line. Therefore, 
the nonuniformity of B r , via the magnetic flux concentration toward 
the magnetic poles (i.e., around the direction of jf), causes an en¬ 
hancement in spindown losses at high inclination angles (a ~ 90°). 

The qualitative behaviour given by eq. (14) has a clear phys¬ 
ical interpretation. From eq. (7) we see that energy flux is mini¬ 
mum near the rotational poles (sin 9 * 0 ) and maximum near the 
rotational equator (sin0 ~ 1). Hence, the spindown losses are min¬ 
imised when the maximum in B, is along the rotational poles and 
does not contribute to the spindown (a = 0). Conversely, the spin- 
down losses are maximised when the maximum in B, lies in the 
rotational equator and contributes to the spindown most efficiently 



@[°] 


Figure 5. The lateral distribution of azimuthally-averaged Z? 2 (which deter¬ 
mines the spindown rate, see eq. 8) due to (i) numerical simulation shown 
with solid curves and due to (ii) the B r distribution for vacuum dipole, 
eq. (10), with dashed curves. Red lines show the results for the aligned and 
solid lines for the orthogonal rotator. Both distributions are normalised to 
have the same magnetic flux. The vacuum dipole approximation accurately 
reproduces the lateral distribution of magnetic field for the orthogonal ro¬ 
tator, as shown with red curves. However, it substantially over-predicts the 
polar field strength and under-predicts the equatorial field strength for the 
aligned rotator. The dotted black line shows an analytic fit ( 18) for the force- 
free aligned rotator simulation results. 


(a = 90°). In other words, the spindown losses are higher than in a 
pure monopole because of angular redistribution of magnetic flux 
from the polar to the equatorial regions. 

From eqs. (10) and (11), the radial magnetic field strength of 
an orthogonal rotator is : 2 

Bf~\r, 6, ip) = Bo(r/ ro ) -2 sin 6 cos ip. (for orthogonal pulsars) 

(15) 

The azimuthally-average of (15) is highly non-uniform laterally 
(see also eq. 13): 

(BX = Q.5Bl(r/r 0 y 4 sin 2 0. (16) 


This nonuniformity implies that the angular distribution of spin- 
down luminosity deviates from the standard expectation, rfL/dw oc 
sin 2 8, eq. (9), and becomes. 



st 2 o 2 

* ^open 

8 ;r 3 cr 2 


sin 4 8, 


(for orthogonal pulsars) 


(17) 


where we used the eq. ( 12 ) to eliminate Bo in favour of <D 0 pen- 


4 APPROXIMATE SOLUTION FOR ALIGNED AND 
ORTHOGONAL FORCE-FREE PULSAR WINDS 

In Section 3.3, we found that magnetic field line bunching around 
the magnetic poles, such as in a vacuum dipole, qualitatively repro¬ 
duces the increased spindown rate of oblique force-free and MHD 
pulsars. This is because most of the spindown emanates near the 
equatorial plane, 8 = n/ 2 , as indicated by the factor sin 3 8 in eqs. ( 8 ) 
and (9). Thus, for a field line to contribute to spindown most effi¬ 
ciently, it needs to be close to the equatorial plane. As we tilt the 


1 Note that since in a dipole the open magnetic field depends on the pulsar 

period, <I> 0 pen x H, the dependence on the period is in 4th power, Tdipole 00 2 This also coincides with B, distribution of the full vacuum solution 

H 4 . (Deutsch 1955). 
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Figure 6. Colour-coded surface distribution of B\ in the numerical force-free solution (red shows high blue low values, see colour bar) for different inclination 
angles, at the surface of a sphere of radius r = 6Rlc- The white contours indicate 91st, 96th, and 99th percentile levels of B r , which was Gaussian-smoothed 
with a standard deviation of 1.5 cells. At all inclination angles the angular distribution of B\ is clearly non-uniform. Even though for an aligned rotator, shown 
in panel (a), the maximum in Bj follows the polar field line (shown with the white asterisk), for other inclination angles (panels b-d) the magnetic peak shifts: 
it moves laterally closer to equatorial plane and shifts in azimuth to a later phase by about </> max » 30° (see also Fig. 7). Note that at nonzero inclination angles 
(panels b-d), the peak in the solution resembles that of the orthogonal rotator (panel d). 


pulsar, its magnetic pole, which is the region of enhanced magnetic 
field, approaches the rotational equator. This leads to an enhance¬ 
ment in the spindown. 

To find out whether this mechanism is at work in the simulated 
pulsar magnetospheres, we compare the angular distributions of B r 
due to a vacuum dipole solution, given by eq. (10), and in a simu¬ 
lation of a force-free magnetosphere. Since we are interested in the 
spindown rate, we will start by looking at an azimuthal average, 
which determines the spindown luminosity (see eq. 8). 

Figure 5 shows the comparison of the vacuum dipole approx¬ 
imation given by eq. (10) and numerical force-free simulations for 
two limiting cases: aligned and orthogonal rotators. The vacuum 
dipole curves are normalised to have the same open magnetic flux 
as in the numerical simulation (at the same obliquity). 

From Fig. 5, we see that the simulation result for the orthog¬ 
onal rotator is well approximated by eq. (16). Thus, the vacuum 
dipole provides an excellent description of an orthogonal rotator 
wind. We adopt eq. (16) as the analytic description of the magneto¬ 
sphere of an orthogonal plasma-filled rotator. 

At the same time. Fig. 5 shows that the vacuum dipole expres¬ 
sion gives a very poor description of the aligned rotator solution: 
whereas in the simulation B, is non-zero in the equatorial region, 
in the vacuum dipole field B, vanishes at the equator. Since the to¬ 
tal magnetic flux is the same, the polar field strength in the vacuum 
model by a factor = 3 exceeds that in the simulation. 

For this reason, we adopt the following fit as the analytic ap¬ 
proximation for the aligned force-free solution: 

B^(8) = S 0 (r/fo) _2 [l + 0.02 sin 6 + 0.22(cos 0-1) 

- 0.07(cos 8 - 1 ) 4 ] x sign cos 8. (18) 

This fit, shown in Fig. 5 with the black dotted line, provides an 
adequate quantitative description of the aligned force-free solution 
(at a representative distance r = 6 Blc by which the solution has 
approximately settled to an asymptotic distribution). We discuss 
the lateral structure of this distribution in Appendix A. 


5 PULSAR WIND STRUCTURE IN MAGNETIC 
COORDINATES 

In §3, we obtained analytic approximations, eqs. (15) and (18), for 
B r (8, (fi) - and thus an approximate description of the pulsar wind - 
for the two limiting cases of aligned and orthogonal rotators. Our 



Figure 7. Positions of the polar field line and maximum in B, are plotted 
with various lines vs the inclination angle, a, as measured at the sphere of 
radius r = 6Rlc- Shown are the polar coordinates of the maximum in mag¬ 
netic field, (0 max , <Pr nax), with error bars indicating the range spanned by the 
outermost white contour in Fig. 6. We also plot the coordinates of the polar 
field line (0 po iar> A>olar), as labeled in the figure. The maximum of the field, 
shown with blue connected triangles, is located near the equatorial plane 
and outruns the polar field line by > 30° for a ~ 30-60°. The polar angle 
of the polar field line, which is shown with the green connected diamonds, 
approximately follows the inclination angle, i.e., the polar field line, despite 
being wrapped around in the azimuthal direction, maintains its polar angle. 
By definition we use the intersection point with the polar field line to set the 
zero point in the azimuth, i.e., we set ^ po i al = 0. 

goal is to obtain the description for pulsar wind at all inclination 
angles. Could we somehow combine these two expressions in or¬ 
der to get a general solution? With this in mind, let us a look at 
the structure of the force-free pulsar winds simulations at various 
inclinations. 

Fig. 6 shows the map of B\ in the numerical force-free solution 
of a pulsar wind on the surface of a sphere of radius ro = 67?lc- By 
this distance, the structure of the wind has frozen out and does not 
change at larger distances. Panels (a)-(d) illustrate solutions for dif¬ 
ferent inclination angles, a = 0,30,60, and 90°, respectively, in two 
ways: (i) via a colour map (see the colour bar) and (ii) with white 
contours that indicate 91, 96, and 99 percentile levels of max |B,.|. 

For each inclination angle, Bj exhibits a well-defined peak. 
The distribution of B r for the aligned rotator, shown in Fig. 6(a), 
peaks at the rotational z-axis, as consistent with Fig. 5. 
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Figure 9. Colour-coded surface distribution of B; in the numerical force-free solution that was averaged out around the magnetic peak. For inclinations 
a > 60°, this axial averaging approach gives a good approximation of solution structure around the peak. However, the averaging washes out the current sheet, 
which is a source of magnetospheric dissipation and potential origin of gamma-ray emission. See Fig. 6 for more details. 



Figure 8. Angular distribution of the (B“) Vm of the numerical force-free 
solution in magnetic coordinates morphs from that of an aligned rotator at 
a = 0° (red line) to that of an orthogonal rotator at a = 90° (blue line). 
The in-between distributions, for intermediate inclination angles, 0° < a < 
90° resemble neither of these two limiting solutions, making it difficult to 
construct an analytic solution for those obliquities. 


How does the situation change for an oblique rotator? We 
show the solution for a 30-degree inclined rotator in Fig. 6(b). The 
white star shows the point at which the polar field line pierces our 
sphere of r = 67?lc- We denote this position as (# po iar, </> po iar) and use 
it to define the zero meridian, and with it our x—z plane. Thus, by 
definition we set = 0. Figs. 6 and 7 show that 9^^ exceeds 
the inclination angle by only a few degrees: this means that as the 
polar field line propagates out and winds up around the rotational 
axis, it does so mostly along a cone of a constant opening angle. 

However, this does not mean that the entire magnetosphere 
tilts rigidly and retains its lateral structure. If it did so, we would 
expect the peak in B, to shift by the inclination angle, i.e. we would 
expect to have 0 max = a. Figures 6(b) and 7 show that this does not 
happen: for a stellar tilt of a = 30°, the peak tilts by three times as 
much, 0 max K 90°. Figs. 6(c) and 7 show that this trend persists for 
higher inclination angles: e.g., for a = 60°, we also find # max « 90°. 
The polar field line and the peak in B r only catch up with each 
other for an orthogonal rotator, a = 0 max = 90°. Therefore, a tilted 
magnetosphere rearranges itself in a non-trivial way, with different 
parts moving relative to each other and causing the maximum of 
radial magnetic field strength to shift toward the equatorial plane 
by more than we would naively expect. 



Figure 10. Angular distribution of the (B x ) v in force-free simulations (solid 
lines) and same numerical simulation result but additionally (C m —averaged 
around the magnetic peak, ((B x )^ m ) (dotted lines) for different obliquities. 
The angle-averaging around the magnetic pole gives a good approximation 
for a = 0° and a > 60° but washes out the non-axisymmetric structure 
around the pole at intermediate inclinations, which manifests itself as the 
smoothing out of the equatorial peak in the dotted orange for 

a = 30°. 


Note that the above picture applies at a rather large distance 
from the stellar surface, r = 67? LC . At much smaller distances, 
e.g., at the surface of the star, the magnetic peak and the polar field 
line coincide. They smoothly diverge from each other with increas¬ 
ing distance from the stellar surface. The solution shown in Fig. 6 
qualitatively persists out to much larger distances; asymptotically 
far away from the star, however, it undergoes small quantitative 
changes. We expect that at larger distances, r > 67?lc. the solution 
does not have much opportunity to rearrange laterally: this is be¬ 
cause the Lorentz factor of the wind scales as y ss R/Rlc » 1, 
and hence we expect that the solution might only be able to change 
on the small angular scales, 59 ~ 1/y ~ Rlc/R :£ 1/6 ~ 10° for 
R > 6A’[,c. 

Despite this magnetospheric rearrangement, the surface dis¬ 
tribution of B;(r 0 ,B, (fi) appears to maintain a rather high degree 
of symmetry around the magnetic peak, (0 max , ^> max ). Are the de¬ 
viations from axisymmetry around the peak important for the an¬ 
gular distribution of emergent luminosity? To test this, we plot in 
Fig. 8 the magnetic field distribution vs the magnetic co-latitude 
9 m (which is counted off from the magnetic peak), averaged over 
the magnetic azimuth i p m (relative to the magnetic peak), (Bj) . 
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We see that as the inclination angle increases, the profile of (B 2 ) 
smoothly transitions from that of an aligned pulsar (red line, a = 
0°) to that of a vacuum dipole, B 2 oc cos 2 8 m (blue line, (B 2 ) ). 

This is another indication that the vacuum dipole magnetic field 
distribution provides a good description of the orthogonal pulsar. 

Fig. 9 shows that after averaging out in tp m , the solution visu¬ 
ally appears very similar to the full force-free simulation. Will this 
approximation give a similar spindown compared to the simula¬ 
tion? To test this, in Fig. Fig. 10 we show the resulting distribution 
of ((B 2 ) ) with dotted lines and the results of the full force-free 

simulation with solid lines. The agreement is excellent at all incli¬ 
nation angles: the differences are generally small and are noticeable 
mostly at small values of a. For instance, at a = 30°, the approxi¬ 
mation smoothes out the otherwise well-defined equatorial peak in 
(B 2 ) v . 

Since most of the spindown emerges near the equatorial plane 
— and these regions are rather well-described by this approxima¬ 
tion — we expect that it provides a decent description of pulsar 
spindown. Flowever, there are two difficulties in using this descrip¬ 
tion in practice. First, it requires knowing the shapes plotted in 
Fig. 8 for all inclination angles. Second, it ignores the presence 
of the equatorial current sheet: in fact, the averaging in <p m “washes 
out” the jump in the magnetic field across the current sheet, which 
is clearly seen as a jump in B, in the bottom right of Fig. 6 (b),(c). 
Since the current sheet could be an important location at which 
dissipation takes place and magnetospheric radio and gamma-ray 
emission is produced, we will try to find a solution that properly 
describes the current sheet. 


6 SEMI-ANALYTICAL SOLUTION TO THE PULSAR 
WIND 

The current sheet structure in the force-free solutions is seen more 
clearly in Fig. 11, which shows the view of the solution down the 
y-axis. The current sheet is clearly seen as the sharp stripe across 
the solution. Apart from a small-scale deformations, the current 
sheet appears to lie in a single plane inclined relative to the x—y 
plane. 

A common description of the magnetospheric current sheet 
is due to Bogovalov (1999) split-monopole model. This approach 
gives a 3D, analytic description of the magnetospheric current sheet 
location under the assumption that the sheet is “painted” on the 
magnetosphere and is passively advected outward at the speed of 
light, as shown in Fig. 12. This model is in rather good agreement 
with the numerical simulation shown in Fig. 11 and captures both 
the inclination and the phase of the current sheet correctly. How¬ 
ever, the split-monopole model assumes a piece-wise constant B r 
outside of the current sheet, and this is clearly not the case in the 
numerical simulations. 

Our goal is to construct an approximate analytic description of 
the 3D pulsar wind that would correctly reproduce both the loca¬ 
tion of the current sheet and the smooth variation of the magnetic 
field everywhere else. In particular, the solution needs to capture 
the gradual decrease of the jump in magnetic field strength across 
the current sheet as the inclination angle increases from a = 30° to 
90°, as seen in Fig. 11 (b)—(d). This might seem like an intractable 
problem. However, as we discussed in §2, the situation is substan¬ 
tially simplified by the fact that to describe the entire wind solution 
it suffices to specify just one function, B, = B°.(8,p), at a suffi¬ 
ciently large radius r 0 » R LC at which the lateral structure of the 
solution “freezes out”. After this, the rest of the solution can be 


unambiguously reconstructed. As in the split-monopole model pul¬ 
sar wind (Bogovalov 1999), this distribution rotates at an angular 
frequency Q of the star in the (/^-direction and “flies out” at essen¬ 
tially the speed of light in the r-direction. The full time-dependent 
solution combines both of these effects: 

B r = B?(a,0, V > + A^)|^J , (19) 

where the phase shift is A tp = £1 • [(r - r 0 )/c - (f - f 0 )], and we 
approximate B e = 0, as the wind is mostly radial at large distances 
(see Figs. 2 and 4). Then, the toroidal magnetic field is given by the 
winding of this radial field eq. (6) (see also Gruzinov 2006). Within 
the approximation = —Eg that we will adopt, the Lorentz factor 
of the flow is 

7 = (g 2 — E 2 ) i/- = §~ = (1 + Q2r2 sin2 ^* /2 ~ Qr sin d ! c 

(Tchekhovskoy et al. 2008). In the simulations, we find that 
B\-El< B 2 « B 2 , E 2 , so the Lorentz factor could be higher 
or lower by a factor ;£ 2 than the estimate (20). For simplicity, we 
will neglect this effect. 

Our goal is to construct an appropriate 2D function, B° r (8,<p), 
that would satisfy the above requirements. As we discussed in §5, 
the structure of the pulsar wind is relatively simple in the two lim¬ 
iting cases: 

(i) aligned rotator (a = 0°). While no analytic solution is 
known, we can write down a simple analytic fitting function, 
eq. (18), for the numerical solution. We will denote this solu¬ 
tion, rotated by an angle a around the v-axis, via B\ l) (r, a, 8, <p) = 
S* 1 ' [r, 8 m (a, 8, tp), <p\, where B^\r,8,p ) is given by eq. (18) and 
8 m (a, 8, ip) by eq (11). 

(ii) orthogonal rotator (a = 90°), whose B, distribution is well- 
described by eq. (15). We will denote this solution as B[~\r, 8, <p). 

Our goal is to describe what happens in between these two 
limiting cases. We established that the position of the current sheet 
is well-described by the split-monopole model, i.e, the current sheet 
lies in the plane that passes through the y-axis and makes an angle 
a, i.e, the inclination angle, with the x—y plane. We can incorporate 
this feature by tilting the solution for the aligned rotator by the 
inclination angle, a, as described in item (i). 

This is clearly not the whole story since this would lead to the 
magnetic peak located at 0 max = a in the x—z plane. However, we 
saw that the position of the magnetic peak is systematically shifted 
toward the equatorial plane, i.e., d max > a. Further, the strength 
of the jump across the current sheet progressively decreases with 
increasing a. Eventually, for a = 90°, the solution becomes the 
orthogonal solution [see item (ii)] without any current sheet. 

On a phenomenological level, this suggests that it would be 
advantageous to mix the tilted solution (i) with the orthogonal one 
(ii). This might work because the magnetic peak of the orthogonal 
rotator is always located at the midplane and will bias the sum to 
peak toward the midplane: 

S“(a, 8, ip) = [w 1 (a)B< 1) (a, 8 , p) + w 2 (a)B (2 \r, 8,<p- </w)] \fr(a), 

( 21 ) 

where if/(a) describes the variation of total open flux vs the incli¬ 
nation angle, as seen in Fig. 3. Here i/> max = 30° is the approx¬ 
imate azimuth of the magnetic peak, as seen in Fig. 7. We can 
determine the weights w i, w 2 , and in the following way. Since 
the orthogonal solution (ii) vanishes at the rotational pole, i.e., 
B^\d = 0, p) = bS- I, (0 = 7T, p) = 0. we choose the value of W 2 
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Figure 11. Colour-coded surface distribution of B\ in the numerical force-free solution at r = 6 Rlc viewed down the y-axis. This is similar to Fig. 6 where the 
same solution is viewed down the magnetic peak. For all inclination angles, except a = 90°, the solution shows a well-defined current sheet whose orientation 
is well-described by the split-monopole model shown in Fig. 12. 



Figure 12. Colour-coded surface distribution of B j: in the split-monopole solution (Bogovalov 1999). The current sheet, in which the radial magnetic field 
vanishes, describes the orientation of the current sheet in the numerical force-free solutions shown in Fig. 6. 


Numerical force-free solution at r = 6 Rlc- 






0.6 
0.4 
0.2 

i|°'° 

Bf. / maxS^ 


Analytical approximation: 



Figure 13. Colour-coded surface distribution of Bj. in the numerical (panels a-d) and analytical (panels e-h) solutions. For convenience of comparison, the top 
row of panels duplicates the numerical solution shown in Fig. 11. The analytical solution, shown in the bottom row of panels, clearly reproduces the main 
features of the numerical solution for all inclination angles. Firstly, the position of the current sheet (the jump in the magnitude of B j:) are in good agreement. 
Secondly, the position and size of the peaks in B y are consistent with each other to within the middle (96th percentile) contour level. Thus, the analytical 
approximation provides an excellent description of the numerical solution. 
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Three-dimensional Pulsar Wind 9 


Numerical force-free solution at r = 6/?lc: 






Analytical approximation: 


0.6 

0.4 

■ 0.2 

|o.o 

Bj. / ma 



Figure 14. Same as in Fig. 13 but centred at the position of the magnetic peak. The analytic model, shown in panels (a)—(d) provides an accurate description 
of the numerical solution at all inclination angles. 


by requiring that we get the polar value of B, from the simula¬ 
tion. Next, we choose Wi such that the weighted sum of the two 
solutions in the square brackets in eq. (21) gives us an inclination- 
independent value of the total magnetic flux. Then, we choose to 
reproduce the dependence of the open magnetic flux on the incli¬ 
nation angle, as seen in the simulations and in Fig. 3. This way we 
find the weights of the aligned and orthogonal components of the 
full solution, respectively, 

Wi(a) = |1 - 2a/n \, (22) 

w 2 (a) = 1 + 0.17 |sin2ar| - wq, (23) 

and the overall normalisation factor that accounts for the variation 
of the open magnetic flux with the inclination angle, 

i //(a) = 1 + 0.2 sin 2 a. (24) 

We have found that eq. (21) with the weights given by eqs. (22)- 
(24) provide a good description of the numerical simulation results. 
We discuss the resulting two-component solution below. 

Figs. 13 and 14 show the comparison of such a two- 
component approximation to the numerical solution from two 
viewing angles: one along the y-axis and the other along the mag¬ 
netic peak, respectively. Fig. 13 shows that our analytic description 
accurately reproduces not only the position of the magnetospheric 
current sheet but also the magnitude of the jump in B 2 across the 
magnetospheric current sheet. This highlights an important differ¬ 
ence from the Bogovalov (1999) split-monopole solution, in which 
B r changes sign across the current sheet but maintains the same 
absolute magnitude, i.e., there is no jump in B 2 across the current 
sheet (see Fig. 12). 

Fig. 14 shows that our analytic approximation is in excel¬ 
lent agreement with the numerical simulation outside of the current 
sheet region. In particular, for all inclination angles we considered, 
the analytic description captures the angular size of the magnetic 
peak and its location (to within 4% in B 2 ). 



Figure 15. Angular distribution of the (S 2 ) v of the numerical force-free 
solution evaluated at r = 6Rlc (solid curves) and our semi-analytical de¬ 
scription (thick dotted curves), for different obliquities: a = 0,30,60, and 
90°, shown with red, orange, green, and blue lines, respectively. Black thin 
dotted line shows the aligned split-monopole solution. The split-monopole 
solution is at the same magnetic flux as the a = 0° aligned force-free so¬ 
lution, shown with the red lines. Clearly, solutions at all inclination angles 
show a substantial degree of non-uniformity in B; that is not present in the 
split-monopole solution. 

Fig. 15 shows that our analytic approximation accurately de¬ 
scribes the azimuthally-averaged radial magnetic field distribution, 
which sets the spindown power according to eq. (8). Importantly, 
it describes the magnetospheric structure not only for aligned and 
orthogonal rotators (by construction), but also for intermediate in¬ 
clination angles. For instance, for a = 30°, it accurately captures 
both the low-latitude hump as well as the high-latitude wings in 
(B 2 )^. Similarly, the solution provides a rather good description 
of the pulsar wind at a = 60°, for which it only slightly over- 
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Figure 16. Angular distribution of the ( dL/dai) v (in arbitrary units) of the 
numerical force-free solution (solid curves) and our semi-analytical de¬ 
scription (thick dotted curves) for different obliquities. Fine dotted black 
curve shows the luminosity due to an aligned split-monopole solution 
and rare dotted black curve shows the luminosity due to an orthogonal 
split-monopole solution. The orthogonal numerical solution has a much 
more pronounced peak, and the peak luminosity exceeds that of the split- 
monopole solution by a factor of 2. This highlights the importance of ac¬ 
counting for the effects of non-uniformity of pulsar wind on the wind lumi¬ 
nosity and its angular distribution. 

predicts (by ~ 5°) the width of the magnetic peak. Note that the 
solutions at all inclination angles are qualitatively different from the 
split-monopole solution, which does not have any non-uniformity 
in 6. For the aligned rotator, this difference emerges from the return 
current present in the magnetosphere: the return current, which is 
not included as part of the split-monopole solution, causes an ad¬ 
ditional j x B force, which is directed away from the midplane, 
and which is compensated by the magnetic pressure gradient force, 
which is caused by a non-uniformity in B, and directed toward the 
midplane (see Appendix A). 

What kind of electromagnetic luminosity does our analytic ap¬ 
proximation imply? How does it compare with the simulated re¬ 
sult? Fig. 16 shows the distribution of luminosity per unit solid an¬ 
gle, (dL/dw/^ oc (B%) sin 2 6 (see eq. 9), as a function of the polar 
angle. The black dotted line shows the expected angular distribu¬ 
tion for a split-monopole model of Bogovalov (1999), dL/da) oc 
sin 2 6. Whereas this dependence is broadly representative of the 
aligned force-free simulation result, which we show with the red 
solid line, it over-predicts the luminosity in the equatorial region 
by as much as 50% at the expense of under-predicting the luminos¬ 
ity at other angles. For other inclination angles, the split-monopole 
solution under-predicts the luminosity by as much as a factor of 2. 
For instance, for a = 90°, the split-monopole solution gives the 
peak value of (dL/da>) ~ 1, much smaller than the simulated peak 
value is s» 2.25. This is despite the open magnetic flux is the same in 
the split-monopole solution as in the numerical simulation. As we 
discussed previously, this excess spindown luminosity in oblique 
rotators is due to the bunching of magnetic field lines toward low- 
latitude regions (see §3). 


7 DISCUSSION AND CONCLUSIONS 

In this work, we analysed the results of 3D, time-dependent numer¬ 
ical simulations of oblique, plasma-filled pulsar magnetospheres 


in the MHD and force-free approximations. Oblique plasma-filled 
pulsars spin down twice as fast as aligned pulsars (Spitkovsky 
2006), i.e., the spindown luminosity L varies with the inclination 
angle a as 

L ~ I + sin - o). (25) 

The reasons for such an increase are not understood. Could this in¬ 
crease be caused by the changes in the amount of open magnetic 
flux? Only to a degree: we found that changes in the open mag¬ 
netic flux account for about 40% of the increase. We argue that the 
other 60% of the increase is due to a nonuniform bunching of mag¬ 
netic field lines toward the equatorial plane. This redistribution al¬ 
lows the field lines to extract the rotational energy more efficiently, 
thereby accounting for the rest of the increase in the spindown lu¬ 
minosity. 

To illustrate this effect, we constructed a toy model, in which 
we approximated the angular distribution of the radial magnetic 
field, B, , by that of a dipole. Dipole magnetic field strength peaks 
at the magnetic poles. As this field distribution rotates with the star, 
the magnetic field lines become wound up around the axis. The 
faster the winding of a field line, the more that field line contributes 
to the spindown. Since the winding is slowest at the magnetic poles 
(there v v = 0) and fastest at the equator {v v is maximum), the 
field lines near the polar magnetic peak do not contribute much 
to the spindown. However, in an inclined rotator, the dipolar peak 
moves closer to the equatorial plane. The field lines near the in¬ 
clined peak contribute to the stellar spindown more efficiently and 
lead to an increase in the spindown of an oblique rotator relative to 
the aligned rotator. We showed that the spindown luminosity in this 
simple model follows the dependence given in eq. (25), which illus¬ 
trates the importance of magnetic field bunching on the spindown 
luminosity. 

The magnetic field of an actual pulsar is not dipolar. Is there 
a quantitative way of describing the non-uniformity of magnetic 
field in the pulsar wind? The good news is that it is sufficient to 
prescribe a single function - a distribution of B, on a sphere of a 
sufficiently large radius - in order to fully recover the asymptotic, 
3D structure of a pulsar wind. In particular, for an aligned rotator, 
the distribution of B, is axisymmetric, and we parametrised it via a 
single scalar function (see eq. 18). 

At first glance, the structure of oblique rotators appears to 
be much more complicated than that of the aligned one. The 
good news is that an orthogonal force-free rotator exhibits a rather 
featureless and simple radial magnetic field distribution, B , oc 
sin 8 cos tp, i.e., the same as of the vacuum dipole that is put on 
its side (see eq. 15)! 

Now, with these two solutions for the asymptotic pulsar wind 
- for the aligned and orthogonal plasma-filled rotators - at hand, 
we combined them into an approximate, two-component solution 
that describes the simulation results at all inclination angles (§6). 

That the full solution is a combination of the two components 
is particularly interesting for the following reason. The magnetic 
field in the aligned component of the solution undergoes a jump 
in the equatorial plane. However, the field remains smooth in the 
orthogonal component. This means that the relative weights with 
which the two solutions enter the full solution tell us about the 
strength of the jump in magnetic field strength and therefore the 
strength of the magneto spheric current sheet. 

As the inclination angle increases, the contribution of the 
aligned solution gradually decreases (eq. 22) and the magneto- 
spheric current sheet weakens, i.e., the jump across the sheet in the 
magnetic field decreases. The current sheet completely disappears 
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for a force-free orthogonal rotator, as seen in Fig. 13(d),(h). This 
suggests that the amount of dissipation in near-orthogonal pulsar 
magnetospheres might be qualitatively smaller than at intermedi¬ 
ate obliquities. This might be the reason for lower dissipation in 
the current sheets found in global PIC simulations (Philippov et al. 
2014). However, the absence of a jump in the magnetic field in a 
force-free solution does not mean that the dissipation is entirely ab¬ 
sent, as found in both recent 3D MHD (Tchekhovskoy et al. 2013) 
and PIC (Philippov et al. 2014) simulations of oblique pulsar mag¬ 
netospheres. 

Whereas force-free and MHD simulations do not include the 
microphysics of the magnetospheric current sheets and reconnec¬ 
tion, they are extremely useful because they can be extended to 
extremely large distances and give us the full, nonlinear structure 
of the pulsar wind electromagnetic field that is otherwise difficult 
to obtain. In particular, they offer us a unique opportunity to de¬ 
termine the smooth structure of the solution outside of the current 
sheets and compute the spindown luminosity, its angular distribu¬ 
tion, and its dependence on the obliquity angle of the pulsar. For 
instance, we find that the current sheets in our analytic description 
the pulsar wind are quite unusual in that the magnitude of |B,.| un¬ 
dergoes a substantial jump across the sheet. Since asymptotically in 
the pulsar wind the total (radial plus toroidal) lab-frame magnetic 
pressure is tightly related to square of the radial magnetic, B 2 , via 
B 2 /87r = (B 2 + B 2 )/8tt ss B 2 ( 1 + Qrr 2 sin 2 6>/c 2 )/87r (see eq. 6), 
the lab-frame magnetic pressure also undergoes a jump across the 
magnetospheric current sheet. 

However, having such a jump is at odds with the expectation 
that the total pressure (which is dominated by magnetic pressure 
outside the sheets), must be the same on both sides of the sheet: oth¬ 
erwise, the current sheet would be pushed to the side with a lower 
pressure until the two pressures came to an equilibrium. There is 
a flaw in this argument: since our simulations are relativistic, we 
need to measure the magnetic pressure in the fluid frame, which is 
different than the pressure in the lab frame because the pulsar wind 
flows outward at a relativistic velocity, with a Lorentz factor y » 1 
relative to the star. Hence, it is the fluid-frame magnetic pressure, 
Pmag = b 2 /8n = B 2 !%ny 2 , that is continuous across the current 
sheet. How can this happen if the lab-frame magnetic pressure is 
discontinuous? This can occur if the fluid elements on the two sides 
of the current sheet move at substantially different Lorentz factors, 
7i ^ 72 , such that even though the lab-frame magnetic pressures 
are different, B 2 /87r t B 2 /8n, the fluid-frame magnetic pressures 
are the same, B 2 /Sny 2 t = b 2 /Sn =» b\!%n = B 2 /&7ty 2 . 

Thus, that the square of the radial magnetic field, B 2 , under¬ 
goes a substantial jump, which we find both in our force-free and 
MHD simulations of oblique pulsar magnetospheres as well as in 
our analytic model, necessarily implies that the fluid elements on 
the two sides of the current sheet move at a relativistic velocity 
relative to each other. This relative motion is independent of the 
internal physics of the current sheet and is a generic feature of 
the large-scale geometry, or topology, of the pulsar magnetosphere. 
This might be interpreted as follows: the reconnection in the cur¬ 
rent sheet is a type of driven, rather than spontaneous, reconnection. 
Thus, in order to accurately predict the properties of current sheet 
reconnection rate, particle acceleration, and emission, local simu¬ 
lations of reconnection need to be constructed for the specific ge¬ 
ometry and type of discontinuity across the current sheet in order 
to properly reflect the driving exerted on the current sheet by the 
global magnetospheric structure (Sironi & Spitkovsky 2014; Guo 
et al. 2014). In particular, local simulations of driven reconnection 


are needed to properly represent the magnetospheric current sheet 
physics. 

Our analytic description of the pulsar magnetosphere provides 
a 3D electromagnetic structure of a pulsar wind near and beyond 
the light cylinder for both aligned and oblique pulsars. This is 
useful for constructing the models of emission from the magneto¬ 
spheric current sheet and modelling pulsar gamma-ray light curves 
(Arka & Dubus 2013; Uzdensky & Spitkovsky 2014). The nonuni¬ 
formity of the pulsar wind and the deviations of lateral structure 
from the split-monopole expectation can be included into and will 
improve the accuracy of the models of pulsar wind nebulae (Ca¬ 
mus et al. 2009; Porth et al. 2014a,b), supernovae remnants, and 
magnetar-powered supernovae (Woosley 2010; Kasen & Bildsten 
2010) and core-collapse gamma-ray bursts (Metzger et al. 2011). 
For all of these systems, if an oblique rotator is a generic situa¬ 
tion, the emerging outflow shows smooth, sinusoidal-like variations 
(e.g., see Fig. 14) separated by jumps at the current sheets. Because 
of these variations, the energy flux is more strongly concentrated 
toward the edges of the emerging outflow than in the standard split- 
monopole approximation for the pulsar wind. 

In the magnetar model of core-collapse gamma-ray bursts, a 
pulsar wind focuses into a pair of tightly collimated jets. If jet 
emission correlates with its power (Nemmen et al. 2012), vari¬ 
ability of the pulsar wind could translate into periodic gamma-ray 
variability that is directly observable. In the particularly clean case 
of an orthogonal rotator, the emerging pulsar wind shows smooth 
sinusoidal-like variation at the rotational period of the neutron star 
that could lead to sinusoidal-like modulation of the high-energy 
emission. That this signature is not routinely seen in the light curves 
of long-duration gamma-ray bursts suggests several possibilities. 
The central compact object could be a black hole or a (nearly) 
aligned rotator and produce no periodic variability. If the central 
object is an oblique rotator, the absence of such variability might 
be suggestive of an emission mechanism that washes out small 
time-scale variability of the central engine. For instance, in sub- 
photospheric gamma-ray emission models (e.g., Giannios 2006, 
2008; Pe’er et al. 2006; Giannios & Spruit 2007; Beloborodov 
2010, 2013), the emission is formed over a wide range of radii, 
which means that the emission is sensitive only to large-scale jet 
structure. How gamma-ray burst emission is produced by inhomo¬ 
geneous jets warrants further investigation. 
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APPENDIX A: ANGULAR DISTRIBUTION OF THE 
POLOIDAL FIELD IN THE ALIGNED PULSAR 
MAGNETOSPHERE 


In this section we show that the angular structure of B r in the 
aligned case can be explained analytically. The force balance 
equation in the axisymmetric steady state can be rewritten as 
(Tchekhovskoy et al. 2008) 


B 2 , - E- 
4nR c 


(R n) 


Bl-E 2 


4 nR 


3 (B 2 -E 2 \ 

dn \ 8tt j 


= 0 , 


(Al) 


where R c is the curvature of poloidal field lines, R is cylindrical 
radius, and n is the unit vector along the local electric field. At 
large distances from the pulsar the poloidal magnetic field is mostly 

radial, so A ~ _ i A an[ j R . n * — cos 6. With this we can rewrite 

’ On r 09 

the force balance as 


b 2 „-e 2 

4 nR r 


■ cos 9- 


4 nR 


l_3_ 

7-de 


B 2 -E 2 
8;r 


= 0 . 


(A2) 


We will evaluate this equation just above the equator, thus neglect¬ 
ing all terms containing cos 9, and R = r sin 6 r. The last term in 
(A2) can be rewritten as 

B 2 - E 2 = 4/ " , + B 2 p (\ - Qrr sin 9 2 /c 2 ), (A3) 

c 2 r 2 sin- 6 p 

where we used the MHD condition E 2 = QrR 2 B 2 /cr, and the Am¬ 
pere’s law, which provides the relation between the toroidal field 
and the enclosed current 7(r, 9) B$ = 2 1/cR. 


1-R 2 /R 2 lc 8 d I_ = 1 - R 2 !R 2 lc OB 2 

p 4nR c c 2 R 2 d9 R 39 ' 

Finally, in the wind region R » R^c we get 

dBl = 8Rl C; dl g2 R 
39 c 2 R 4 d9 "4 tiR c ' 


(A4) 


(A5) 


We find that the first term dominates at R > 27? ic and is nega¬ 
tive just above the equator because of the presence of the volume- 

Q g2 

distributed return current. Finally, we obtain < 0 as it is in the 
simulation. 


APPENDIX B: FINITE MAGNETISATION EFFECTS 

Note that in many of the systems, such as in PWNe and SNRs, 
pulsar wind propagates for many decades in radius, so that eventu¬ 
ally the flow becomes super-fast magnetosonic. This happens at the 
fast magnetosonic surface, at which the flow reaches Lorentz factor 
7 f ~ 7 max> where y max is the terminal Lorentz factor that the outflow 
would achieve if all of its electromagnetic energy were converted 
into kinetic energy flux (Beskin et al. 1998; Tchekhovskoy et al. 
2009). For pulsars, theoretical estimates suggest y max ~ 10 4 — 10 6 
(e.g., Kennel & Coroniti 1984a,b). Since in a force-free flow the 
Lorentz factor increases approximately linearly with distance from 
the star, y » r sin 9/R^c, it crosses the fast surface at r F ~ ~ 

(20-100)T? LC . Beyond this distance, y ~ /r 1 / 3 /u 1 / 3 (r/r F ) (see Ap¬ 
pendix B in Tchekhovskoy et al. 2009). The fast surface is located 
at a very short distance compared to those of interest in PWNe and 
SNRs, which means that finite-magnetisation, MHD effects will set 
in long before the pulsar wind starts to radiate due to the interaction 
with the ambient medium in these astrophysical systems. Will this 
lead to large deviations away from the force-free structure that we 
have considered? We are interested in the field lines that make a 


small angle with the midplane, 9' = 7 r /2 - 9, since it is these field 
lines that carry most of the energy flux. Fortunately, once the out¬ 
flow accelerates up to a large Lorentz factor, it is extremely difficult 
to reorient it. In fact, the relative change in the direction of the flow 
A ff Iff = y/y max ~ jU ' 2/3 ln 1 / 3 (r/rp) <K 1 , even for a really small 
H ~ 10 2 and a really large r = 10 1( V F . The only place where the 
finite magnetisation MHD effects might make a difference are the 
polar regions of the flow, which, however, are energetically unim¬ 
portant. 
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